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Two-photon photoassociation of hot magnesium atoms by femtosecond laser pulses is studied from first prin- 
ciples, combining ab initio quantum chemistry and molecular quantum dynamics. This theoretical framework 
allows for rationalizing the generation of molecular rovibrational coherence from thermally hot atoms [L. Ry- 
bak et al, Phys. Rev. Lett. 107 , 273001 (2011). larXiv:1107.4755] . The coupled cluster and multi-reference 
configuration interaction frameworks are used to calculate the relevant potential energy curves, one-photon 
and two-photon transition matrix elements, dynamical Stark shifts, as well as spin-orbit couplings and non- 
adiabatic radial couplings. Random phase thermal wavefunctions are employed to model the thermal ensemble 
of hot colliding atoms. Comparing three different choices of random phase thermal wavefunctions, the free 
propagation approach is found to have the fastest initial convergence for the photoassociation yield. When 
further refinement is required the eigenvalue approach is superior. The interaction of the colliding atoms with 
two time-delayed femtosecond laser pulses is modeled non-perturbatively to account for strong-field effects 
observed in the experiment. Good agreement between experimental and theoretical results is obtained. 



I. INTRODUCTION 

Molecules can be assembled from atoms using laser 
light. This process is termed photoassociation. With 
the advent of femtosecond lasers and pulse shaping tech- 
niques, photoassociation became a natural candidate for 
coherent control of a binary reaction. Coherent control 
had been conceived as a method to determine the fate of 
chemical reactions using laser fields X The basic idea is 
to employ interference of matter waves to constructively 
enhance a desired outcome while destructively suppress- 
ing all undesired alternatives^ Control is exerted by 
shaping the laser pulses, the simplest control knobs be- 
ing time delays and phase differences^ Over the last two 
decades, the field of coherent control has developed sig- 
nificantly both theoretically and experimentally^— How- 
ever, a critical examination of the achievements reveals 
that successful control has been demonstrated almost ex- 
clusively for unimolecular processes such as ionization, 
dissociation and fragmentation. It is natural to ask why 
the reverse process of controlling binary reactions^— is 
so much more difficult. 

The main difference between unimolecular processes 
and a binary reaction lies in the initial state - a single or 
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few well-defined bound quantum states vs an incoherent 
continuum of scattering states^ For a binary reaction, 
the nature of the scattering continuum is mainly deter- 
mined by the temperature of the reactants. As temper- 
ature decreases, higher partial waves arc frozen out. At 
the very low temperatures of ultracold gases, the scatter- 
ing energy of atom pairs is so low that the rotational bar- 
rier cannot be passed, and the scattering becomes purely 
s-wavei 2 ^ In this regime, the reactants are pre-correlated 
due to quantum threshold effects 2 ^ and the effect of scat- 
tering resonances is particularly pronounced i^2r— At a 
temperature of about 100 /iK, photoassociation with fem- 
tosecond laser pulses has been demonstrated^ Coherent 
transient Rabi oscillations were observed as the promi- 
nent feature in the pump-probe spectra. The transients 
are due to long tails of the pulses caused by a sharp 
spectral cut which is necessary to avoid excitation into 
unbound states. 25,26 This pinpoints to the fact that the 
large spectral bandwidth of a femtosecond pulse is not 
well adapted to one-photon photoassociation at ultralow 
temperatures which requires narrow-band excitation*^ 

The situation changes completely for high tempera- 
tures where the scattering states can penetrate rotational 
barriers due to the large translational kinetic energy. The 
association process is then likely to happen at short in- 
ternuclear distance close to the inner turning point and 
for highly excited rotational states. In this case, the 
large spectral bandwidth of femtosecond laser pulses is 
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ideally adapted to both the broad thermal width of the 
ensemble of scattering states and the depth of the elec- 
tronically excited state potential in which molecules are 
formed. The disadvantage of this setting is that the ini- 
tial state is completely incoherent, impeding control of 
the photoreaction. Photoassociation with femtosecond 
laser pulses was first demonstrated under these condi- 
tions, employing a one-photon transition in the UVJ^ 
Subsequent to the photoassociation, coherent rotational 
motion of the molecules was observed*^ We have recently 
demonstrated generation of both rotational and vibra- 
tional coherences by two-photon femtosecond photoasso- 
ciation of hot atomsi 28 ' 29 This is a crucial step toward the 
coherent control of photoinduced binary reactions since 
the fate of bond making and breaking is determined by 
the vibrational motion. 

Employing multi-photon transitions comes with sev- 
eral advantages: The class of molecules that can be pho- 
toassociated by near-IR/ visible femtosecond laser pulses 
is significantly larger for multi-photon transitions com- 
pared to one-photon excitation. Femtosecond laser tech- 
nology is most advanced in the near-IR spectral re- 
gion. Due to the different selection rules, different elec- 
tronic states become accessible for multi-photon transi- 
tions compared to one-photon excitation. Control strate- 
gies differ for multi-photon and one-photon excitation. In 
particular, large dynamic Stark shifts and an extended 
manifold of quantum pathways that can be interfered 
come into play for multi-photon excitation^ The theo- 
retical description needs to account for these strong-field 
effects. 

We have constructed a comprehensive theoretical 
model from first principles to describe the experiment in 
which magnesium atoms in a heated cell arc photoasso- 
ciated by femtosecond laser pulsesi 28 ' 29 It is summarized 
in Figure [TJ Magnesium in its electronic ground state is a 
closed shell atom. Its ground electronic potential, X lv J+, 
therefore displays only a weak van der Waals attractive 
well. A femtosecond pulse of 100 fs transform-limited du- 
ration with a central wavelength of ~ 840 nm promotes 
an electron to the 7r orbital. This two-photon transition is 
driven since a wavelength of 840 nm is far from any one- 
photon resonance both for magnesium atoms and Mg2 
molecules, cf. Fig. [TJ Upon excitation, a strong chem- 
ical bond is formed in the (l) 1 II g state with a binding 
energy of ~ 1.8 cV or, equivalently, ~ 14500 cm -1 . A 
time-delayed femtosecond pulse probes the excited Mg2 
molecule by inducing a one-photon transition to a higher 
excited electronic state ( 1 I1 U ). This state has a strong 
one-photon transition back to the ground state. The 
corresponding experimental observable is the intensity 
of the resulting UV fluorescence (~ 290 nm), measured 
as a function of the pump-probe time delay. An oscil- 
lating signal is a manifestation of coherent rovibrational 
dynamics in the 1 II S statei 28 i 29 

The correct description of the thermal initial state is 
crucial to capture the generation of coherence out of an 
incoherent ensemble. The density operator, p T , describ- 
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FIG. 1. Potential energy curves of the electronic states in- 
volved in the two-photon photoassociation probed by a time- 
delayed pulse. The shaded region indicates the vibrational 
band populated after photoassociation. 



ing the initial state of hot atom pairs at temperature, T, 
is constructed by a thermal average over suitable basis 
functions. Since no dissipative processes occur on the 
sub-picosecond timcscale of the experiment, the coher- 
ent time evolution of the density operator is most easily 
carried out by propagating the basis functions. Expec- 
tation values are obtained by thermally averaging the 
corresponding operator over the propagated basis func- 
tions. A numerically efficient description of the initial 
thermal ensemble is thus essential to facilitate the time- 
dependent simulations. This is the main subject of our 
present work. 

The paper is organized as follows. Section [TTJ iden- 
tifies the relevant electronic states and presents their 
potential energy curves, transition matrix elements and 
non-adiabatic couplings. These are obtained employing 
highly accurate state of the art ab initio methods. An 
effective description of the thermal ensemble of transla- 
tionally and rotationally hot atom pairs in their elec- 
tronic ground state based on random phase thermal 
wavepackets is developed in Section IIIII Three different 
choices of basis functions are discussed. Section [IV] intro- 
duces the light-matter interaction. The time-dependent 
Schrodinger equation is solved separately for the pump 
and probe pulses. Convergence of the photoassociation 
probability is studied in Section [V] for the different ther- 
mal averaging procedures, and the role of shape reso- 
nances is discussed. Section [VTJ investigates the gener- 
ation of coherence in terms of the quantum purity and 
a dynamical coherence measure. Finally, we conclude in 
Section ELD 
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II. AB INITIO POTENTIAL ENERGY CURVES, 
TRANSITION DIPOLE MOMENTS, AND 
NONADIABATIC COUPLING MATRIX ELEMENTS 

State-of-the-art ab initio techniques have been applied 
to compute the potential energy curves of the magne- 
sium dimcr in the Born-Oppcnhcimcr approximation. In 
all calculations the aug-cc-pVQZ basis set of quadruple 
zeta quality as the atomic basis for Mg were used. In ad- 
dition, this basis set was augmented by the set of bond 
functions consisting of [3s3p2d2flglh] functions placed 
in the middle of the Mg2 dimer bond. All potential en- 
ergy curves were obtained by a supermolecule method 
and the Boys and Bcrnardi scheme was used to correct 
for the basis-set superposition error.— 

The ground X 1 !!^ state potential was computed with 
the coupled cluster method restricted to single, double, 
and noniterative triple excitations, CCSD(T). For the ex- 
cited 1 II g and (l) 1 !^ states we used the linear response 
theory (equation of motion approach) within the coupled- 
cluster singles and doubles framework, LRCCSD. The 
potential energy curve of the excited (2) 1 II ll state in the 
region of the minimum of the potential was also obtained 
with the LRCCSD method, whereas at larger internu- 
clear distances this potential energy curve was repre- 
sented by the multipole expansion with electrostatic and 
dispersion terms C n /R n up to and including n = 10. The 
long-range coefficients C„ were obtained within the mul- 
tireferencc configuration interaction method restricted to 
single and double excitations, MRCI, with a large ac- 
tive space. The latter procedure was necessary since the 
(2) 1 n u state dissociates into Mg( 3 P)+Mg( 3 P) atoms and 
cannot be asymptotically described by a single Slater de- 
terminant. The energy of the separated atoms was set 
equal to the experimental value for each electronic state, 
although the atomic excitation energies obtained from 
the LRCCSD calculations were very accurate and for the 
lowest X P state the deviation from the experimental val- 
ues was approximately 100 cm -1 . 

A high accuracy of the computed potential energy 
curves is confirmed by an excellent agreement of the 
theoretical dissociation energy for the ground A^E^" 
state (Do =403. lem -1 ) with the experimental value 
(Do =404.1±0.5cm _1 )i22 Moreover, the number of 
bound vibrational states for J = supported by the elec- 
tronic ground state agrees with the experimental num- 
ber N u = 19. Spectroscopic parameters of another 
experimentally observed state, ^4 1 S+, also agree with 
our values, for the well position within 0.07 Bohr, while 
the binding energy (D e =9427cm™ 1 ) is only 0.4% higher 
than the experimental value (D e =9387cm _1 )< 3 ^ Also the 
root mean square deviation of the rovibrational levels 
computed with the potential energy curves from the 
CCSD(T) and LRCCSD calculations for the ground and 
A states were 1.3 cm -1 and 30 cm -1 , respectively, i.e. 
0.3% of the potential well depth. Such a good agreement 
of the calculations with the available experimental data 
strongly suggests that we can expect a comparable level 



of accuracy for other computed potential energy curves 
and molecular properties. More details on the ab initio 
calculations and spectroscopic characteristics of various 
electronic states of Mg2 will be reported elsewhere^ 3 - 

The electric transition dipole moments between states 
i and / \ii = where the electric dipole oper- 

ator /ij = 7*j is given by the ith component of the posi- 
tion vector and ^i/f are the wave functions for the ini- 
tial and final states, respectively, were computed as the 
first residue of the LRCCSD linear response function for 
A 1 !]^ - , 1 LT S , and 1 II tl (l) states, whereas for transitions 
to (2) 1 LI„ state MRCI method was applied. 

Nonadiabatic radial and angular coupling matrix ele- 
ments, as well as the spin-orbit coupling matrix elements 
have been evaluated with the MRCI method. The adia- 
batic (1) 1 LI U and (2) 1 H U states are strongly coupled by 
the nuclear momentum operator T# = iV k/a 4 - There- 
fore to include this non-adiabatic coupling in the model 
one has to switch to the diabatic representation, i.e. to 
using diabatic potential energy curves and coupling po- 
tentials: 



V(R) = 



_(v 1 \{R) v 12 (R) 



V 21 (R) V 2 d 2 (R) 



(1) 



where Vj^ and V 2 % are the diagonal diabatic potentials, 
while V\ 2 (R) = V 2 i(R) is the coupling term. The poten- 
tial curves corresponding to the diabatic potentials 
and V22 ma y cross, and are obtained from the following 
expression, 

vro-A-'w Vm ° (r) )mh), 

(2) 

where V( 1 / 2 ) 1 n„ (R) are the adiabatic potential energy 
curves for the (l) 1 !^ and (2) 1 U U states, and the matrix 
A(R) takes the form 



MR) 



cos C(R) sin ((R) 
-sinC(-R) cos COR) 



where 



C(R) 



T(R')dR' . 



The nonadiabatic radial coupling t(R) is given by: 

d 



t(R) = { * ( i)i n „ 



dR 



(2)in t 



(3) 



(4) 



(5) 



Diabatic molecular wave functions are constructed 
from the adiabatic ones by applying the transformation 



*d2 



MR) 



*(2) I n„ 



(6) 



so the corresponding transition moments form the state 
i to the diabatic states can be obtained from those in the 
adiabatic basis 



= MR) 



(^i\X\^ {2} 2 nu 



(7) 
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Weak spin-orbit coupling between the (l) 1 !^, ( 1 ) 3 II S , 
and (l) 3 £ g states was also taken into account in our dy- 
namical model. The spin-orbit matrix elements relevant 
for our work read: 

Wi(iZ) = (* ( i)in 8 |irso|*(i)3 Ee ) J (8) 
W 2 (R) = (* (lV n 9 \H so \y (ir n 9 ), (9) 



W 3 (R) = (*(3)inJffso|*(i)3s 9 >, 



(10) 



where i?so the spin-orbit coupling Hamiltonian in the 
Breit-Pauli approximation including all one- and two- 
electron terms. The spin-orbit coupling matrix elements 
were computed with the multireference configuration in- 
teraction method restricted to single and double excita- 
tions (MRCI). 

The two-photon electric dipole coupling, 

X (uL,t,R) = \\S{t)\ 2 Y,E l E J M[^{u JL ,R) 1 (11) 



is given by the tensor elements of the two-photon elec- 
tric transition dipole moment between the ground and 
excited / state, 



III. QUANTUM DYNAMICAL DESCRIPTION OF A 
THERMAL ENSEMBLE 

The initial state for photoassociation is given by the 
ensemble of magnesium atoms in the heated cell which 
interact via the X X T,^ electronic ground state potential. 
Assuming equilibrium, the initial state is represented by 
the canonical density operator for N atoms held in a vol- 
ume V at temperature T. Due to the moderate density 
in a heat pipe, the description can be restricted to atom 
pairs. The density operator for N/2 atom pairs is then 
obtained from that for a single atom pair,— p T (t = 0), 
which is expanded into a suitable complete orthonormal 
basis. Thermally averaged time-dependent expectation 
values of an observable A are calculated according to 



M 



n#0 



(/l/y™)H£j|0) , </|Aj|n><n|Ai|0> 



(12) 

where lul — hcj\h- In practice, the two-photon tran- 
sition moment given by the expression above can be 
obtained as a residue of the cubic response function.— 
For transitions between the X 1 ^ and states, 

M(^°{ujl,R) was computed as a residue of the coupled 
cluster cubic response function with electric dipole opera- 
tors and wave functions within the CCSD frameworkj 35 ' 36 
The Stark shift of k state can be expressed in terms of 
the elements of the dynamic electric dipole polarizability 
tensor, 



4 K, t, R) = - - 1 S(t) | 2 £ EiEj a% K , R) , 



(13) 



where the tensor elements of the dynamic polarizability 
are given by 



a^(oj L ,R) 



£ 



(k\fk\n)(n\jlj\k) (k\fij\n}(n\p, t \k) 



(14) 

The tensor elements of the electric dipole polarizability 
of the ground X 1 ^ state were obtained as the coupled 
cluster linear response function with electric dipole op- 
erators and wave functions within CCSD framework^ 
Dynamic polarizabilitics of the exited states were com- 
puted as double residues of the coupled cluster cubic re- 
sponse function with electric dipole operators and wave 
functions within CCSD framework i 38 ' 39 

The CCSD(T) and CCSD calculations (including re- 
sponse functions calculations) were performed with the 
DALTON program,— while MRCI calculations were done 
with the molpro suite of codes4i 



(A(t)) T = Tr[Ap T (t)] 



(15) 



with the time evolution of p T (t) given by p T (t) = 
0(t,0)p T (t = 0)0 (t, 0) starting from the initial state 

Pr(t = 0) = \e~^ , 

where f3 = l/fc^T and Z = Tr[e _/3H ] the partition func- 
tion (since we are not interested in calculating absolute 
numbers, normalization by Z will be omitted below). For 
a thermal, i.e. incoherent, initial state, undergoing co- 
herent time evolution, the dynamics can be captured by 
solving the Schrodinger equation for each basis function. 
Solving the Liouville von-Neumann equation for the den- 
sity operator is not necessary if dissipation can be ne- 
glected on the timescales considered. Thermally aver- 
aged expectation values are then calculated by properly 
summing over the expectation values obtained from the 
propageted basis states^ 

Since many scattering states in a broad distribution 
of rotational quantum numbers are thermally populated, 
the approach of propagating all thermally populated ba- 
sis states directly^ becomes numerically expensive. In- 
stead, an effective description of the thermal ensemble 
of scattering atoms is obtained by implementing ther- 
mal averaging as an average over realizations of random 
phases. It makes use of thermal random wavefunctions, 
\ipnj}- Here, the index k labels a set of random phases 
and jn the basis states. Choosing an arbitrary complete 
orthonormal basis set, and given that 



i N 



fe=i 



for random phases 0*, 0g and N large, an expansion into 



unity^ 3 -, 



random phase wavefunctions yields a representation of 

N N 

i = ^Ei**><**i = ^EE e<(BS " 9jl) i«>^i 

fe=l a/3 
1 N 



k=l 



fc=l a/3 
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where = e i9 -\a} and \$> k ) = £ a e i9 »|<*)- Iden - 

tifying a set of suitable basis states \a), only rovibra- 
tional states of the electronic ground state Hamiltonian 
H g need be considered since no electronic excitations are 
excited thermally. A separation of rotational and vibra- 
tional dynamics and subsequent expansion into partial 
waves is natural. Neglecting ro-vibrational couplings, 
the Hamiltonian for the vibrational motion is represented 
on an equidistant grid (Nr = 1024, Rmin = 3.0 ao, 
Rmax = 40.0 ao) for each partial wave J. In the fol- 
lowing, we discuss three possible basis sets for the vibra- 
tional Hamiltonian, from which random phase functions, 
are generated. All three possibilities will lead when av- 
eraged to the initial thermal ensemble corresponding to 
the experiment. While formally equivalent, convergence 
of the thermal averages with respect to the number of 
required basis functions differs significantly for the three 
representations . 



A. Grid-based random phase approach 

The simplest but, as it turns out, most inefficient ap- 
proach uses, for each partial wave J, the coordinate ba- 
sis of ^-functions localized at each grid point, R, 1 = 
J2 R \R, J)(R, J\M A random phase wavefunction is ob- 
tained by multiplying each basis state with a different 
random phase, 9rj, 



I* 



R.J, 



= e «. J 



R,J), 



where k labels one realization of random phases, {Or^j}. 
The resulting wavefunction, (R\^ R j), has constant am- 
plitude and a different random phase at each R. In this 
basis, the initial density operator is obtained by prop- 
agating each \^ k R J ) under H g 7 = f + V g (R) -J- J(J+1) 



in imaginary time, r = i^f3. This yields 



r 2mff 
k 



R,J/T 



\ty R j) and thus the initial density operator, 



N 



Pr(t = 0) = \e-^e-^ 1 £ £ (2 J + 1)(2J' + l)^'"^''^, J){R >, f\ 



1 1 

ZN 



k=l R,J,R'J' 



N 



£e-* H .|n, J K*&>-* n ' 



k=l 



r 



A Chebychev propagator is employed for the imaginary states for the dynamics under the probe pulse. Thermally 

time propagation^ to obtain the N thermal random averaged time-dependent expectation values are obtained 

phase wavefunctions \^ k R j)t which represent the initial from Eq. (|15[) . using cyclic permutations under the trace, 

I 



Tr 



1 



N 



Ap T (*) = T7 E E( 2 J + l) 2 (*l,/l AU(i, 0)p T (t = 0)0 + (t, 0)|* 



N 
1 1 



k \ 

R,j) 



fe=l R,J 
N 



v E E( 2J + l) 2 (*l,7|e- f fl « + (i, 0)A0(t, 0)e"K \* k R>J ) 

k=l R. J 

N 



Jf E E( 2J + T(^j(t)\A\^j(t)) T , 



(17) 



fe=l R,J 



where \^f R j{t))T = U(f, 0)\^ R j)t are the propagated 
random-phase thermal wavefunctions. Note that while 
I*™ j)t has zero components on all electronic states ex- 
cept the ground state, \^ R j(t > 0)}t will be non-zero for 
all electronic states due to the interaction with the field. 
Typical expectation values are the excited state popu- 
lation after the pump pulse, possibly J-resolved. The 
corresponding operators are the projectors onto the elec- 



tronically excited state, P e = |e)(e| and P e , i.e., 



Z(P e J (tf)) = E E( 2J + l) 2 Ke|*^j(*/)>T| 2 ■ 

fe=l R 

The convergence of this approach is slow because there 
is no preselection of those basis states that are most rel- 
evant in the thermal ensemble. 
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B. Eigenfunction-based random phase approach 

A preselection of the relevant states becomes possi- 
ble by choosing the eigenbasis of H g and evaluating the 
trace only for basis states with sufficiently large thermal 
weights, e 2 9 > e. The cigcnfunction-bascd random- 
phase wavcfunctions are given by 



where J denotes the partial waves and n is the vibra- 
tional quantum number in the bound part of the spec- 
trum of H g or, respectively, the label of box-discretized 
continuum states. It is straightforward to evaluate the 
representation of the initial density operator in this ba- 
sis, 



\n,J) 



p T (t = 0) = — e"2 H "e" 



7^E E (2J+1)(2J' + 1) 



A' 



N 



^E E (2J+l)(2J' + l)e^.- e -.-)|n,J)(n',J'| 



(18) 



k=l n,J,n',J' 



Z N 



k—1 n,J,n' \J' 



\n,J){n',J'\ 



r 



where E n ^j are the eigenvalues of the partial wave 

ground state Hamiltonian, fi g . Thermally averaged 
time-dependent expectation values are calculated anal- 
ogously to Eq. ([TT]) . 



Tr 



Ap T (i) 



1 

N 



N 



^^(2J+l) 2 T « J (0|A|^ J (t)) T 



(19) 



k=l n,J 



r 



where the time-dependent wavcfunctions are the (real- 
time) propagated, Boltzmann-weightcd cigenstates with 
random phases, 9^ j, 



l*5,j(*)>T = 0(t,0)|** r)T = e 



'0(t,0)\n,J) 



The eigenfunction-based random phase approach re- 
quires diagonalization of the partial wave ground state 

Hamiltonians, H g . Depending on the time required for 
the propagation of each basis state, this effort may very 
well be paid off by the much smaller number of basis 
states that need to be propagated. 



C. Gaussian random phase wave packets / freely 
propagated random phase wavepackets 

The third approach avoids diagonalization of the par- 
tial wave ground state Hamiltonians, H g , approximating 
them by the kinetic energy, T, only. This approximation 
is valid at high temperatures where the kinetic energy of 
scattering atoms is much larger than their potential en- 
ergy due to the interaction. It employs free thermal wave- 
functions which are preselected according to their Boltz- 



mann weight^ For each partial wave J, free thermal 
wavcfunctions are obtained from (5-functions in momen- 
tum space, 1 = ^2 P \P 7 J)(P, J\ 7 with a random phases, 
8p,Ji 



I* 



P,J) = e 



■ J \P,J) 



Applying e 2 "» e Wp - J \P, J) « e~H^e w '^ J \P, J) and 
Fourier transforming to coordinate space, thermal Gaus- 



I* 



(«-Bor 



\R,J), are obtained. Their 

width is given by (Jr,t = ^/<^p,t = y2m//S = yJirnksT , 
i.e., it matches the Boltzmann distribution. The random 
phase, #p,j, translates into the position of the Gaussian, 
Rq. The density operator can thus be constructed from 
Gaussian functions in coordinate space by either posi- 
tioning them at various points Rq and averaging over Rq, 
or propagating one Gaussian function with a single i?o 

under H g for randomly chosen times, t^ rop and averaging 
over those times. In the latter procedure, the probabil- 
ity amplitude of the two atoms is confined by the repul- 
sive wall of the potential at shorter internuclear distances 
and an artificial reflecting barrier at long distance. The 
wavefunction then keeps on oscillating between the inner 
repulsive wall of the potential and the outer reflecting 
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barrier. The convergence of this method with respect to 
the photassocaion yield is initially very fast, only a few 
realizations are sufficient. The drawback of the proce- 
dure is that only the free part of the initial wavefunctions 

I 



is represented leaving out the interaction energy of the 
true scattering states. Averaging over propagation times, 



thermal expectation values are obtained by 



Tr 



(20) 



fc=i 



r 



where 



'(t)h 



Q(t,0)\Hf t r v ) T and |tf t /"»>T 



\R, J). 



IV. INTERACTION OF THE DIATOMS WITH PUMP 
AND PROBE PULSES 



and thermally averaging the solutions. Absorption of the 
pump and probe pulses is considered separately neglect- 
ing coherent effects when pump and probe overlap in time 
which arc unlikely to play any role here. The Hamilto- 
nian describing the interaction of the diatom with the 
pump laser pulse is given by 



The interaction of the atom pair with the laser fields 
of the pump and probe pulses is simulated by solving 
time-dependent Schrodinger equations, 



it) 



dt 



A(*)l<j(*)>. 



(21) 



J 



Hi 



M 



( H X1E + x*{t,R)e-^ \ 

X (t,R)e^W H' ( 7 1)lno Wi{R) W 2 {R) 




V 







Wi{R) 
W 2 (R) 



Ha)3 Sg W 3 (R) 
W 3 (R) H J (irn J 



(22) 



where H. 



J a(e) 



(t,R) is the nuclear Hamiltonian of elec- 



tronic state a, 



H' 7 



V a (R) 



J(J + 1) 



2mR 



2 



-2 



with T = P /2m the vibrational kinetic energy, m the 
reduced mass, V a (R) the potential energy curve of elec- 
tronic state a, and J the rotational quantum number. 
The two-photon coupling between the A 1 E+ (g) and 
(l^Hg (e) states is denoted by %(t, ii). The strong 

I 



laser field driving the two-photon transitions may lead 
to non- negligible dynamic Stark shifts in the A 1 !]^ (g) 
and (l) 1 !^ (e) states. Wi(R) are the spin-orbit coupling 
matrix elements, cf. Eqs. (|5|-(fTr 



The probe pulse promotes population from the (l^Tlg 
excited state to the non-adiabatically coupled (l) 1 !^ and 
(2) 1 n tl states. Taking again the spin-orbit coupling of 
(l) 1 n s excited state to the (1) 3 £ 9 and (l) 3 n g states into 
account, the Hamiltonian describing the interaction of 
the photoassociated molecules with the probe pulse reads 



HprobeM 



E pr (t)[n(R) 
E pr (t)^ 2 (R) 
Wi{R) 
W 2 {R) 



E pr {t)m{R) E pr (t)^ 2 (R) Wx{R) W 2 (R) 

Ai V 12 (R) 

V 12 (R) H 22 









J 

(i) 3 s 9 





W 3 (R) 



(23) 



W 3 (R) H (1)3ng y 



where Hi/ 2 (R) are the transition dipole moments between 



the (l) 1 !!^ state and the l IL u states. 



T + V£(R) 



J{j+i) 

2mK A ■ 



1 or 2, where V$(R) is the diagonal diabatic 
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potential, cf. Eq. (TJJ. Vi2(R) is the coupling between the 
and 1 II„(2) states in the diabatic representation, 
and E pr (t) denotes the electric held of the probe laser 
pulse. 

The photoassociation is studied numerically by solv- 
ing the time-dependent Schrodinger equations, Eq. (|2ip . 



for the Hamiltonian H pump , Eq. (|22j) . with a Chcbychev 
propagator^ The pump laser pulse excites a small por- 
tion of the incoherent ensemble of ground state atom 
pairs. This action corresponds to a filtering and might 

I 



thus lead to higher purity and coherence of the photoas- 
sociated molecules. In order to study the excited state 
purity, 



V e (t) = Tr[p 2 T Jt)} 



(24) 



and the dynamical coherence measure, C e (t), we con- 
struct the normalized density operator of electronic state 

e > PT,e(*)i 



N 



PT,e(t) = 7p^T^E E (2J + l)(2J / + l)P e |*S, J (t))TT<** V /(t)|P e 

\ ey 'I k=ln,J,n'J' 



In the grid representation using Nr grid points, this becomes a matrix of size Nr x Nr, 

1 1 N 



k—l n,J,n' , J 7 



where ^ e nJ (R,t) = (R,e\^ tJ (t)) T is the excited state 
projection of the propagated fcth thermal random phase 
wavefunction. Since we expect to populate only a limited 
number of excited state eigenfunctions, say N m , it is com- 
putationally advantageous to transform the excited state 
component of the propagated thermal wavefunctions into 



the eigenbasis of the electronically excited state, 



with 



(t) 



^j(R,t)^ m (R)dR. 



The resulting density matrix, 



PT,e (*) = 



l l 



'PJt)) N ^ 

' ey 11 fc=l n,J,n'J 



E (2J+l)(2J' + l)c: 



k,n,J 



(t)c k S' J "(t), 



is only of size N m x N m and can more efficiently be 
squared. Moreover, it lends itself naturally to the evalu- 
ation of the dynamical coherence measure. In the eigen- 
basis, we can easily decompose the density operator into 
its static (diagonal) and dynamic (off-diagonal) part, 
P = Pstat + Pdyn- Such a decomposition has been moti- 
vated in the study of dissipative processes, in particular 
by the fact that pure dephasing does not alter the static 



part 



47.- 



The dynamical coherence measure, 



C e (t) = Tr[p 2 T . eAyn (t)} , 



(25) 



V. CONVERGENCE OF THE THERMAL AVERAGING 
PROCEDURES: PHOTOASSOCIATION PROBABILITY 



The convergence with respect to the number of random 
phase realizations, N, is checked with respect to Boltz- 
mann weighting as follows. Each random phase wave 

t k 

function, \^\j)t, |*„,j)t, or \^J rop ) T , is projected 
onto the eigenstates, | n, J) , of the electronic ground state. 
The resulting overlaps are averaged over all possible re- 
alizations, k, 



captures the part of the purity that arises from the dy- 
namical part of the density operator ! 47 ' 48 



fe=i 



1 



JV 



(26) 
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FIG. 2. Logarithm of the overlaps C^j, cf. Eq. (HHJt, 

vs en- 
ergy, E g n j (J=0) for the grid-based random phase approach. 
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FIG. 3. Initial thermal density PT, g {R) / R 2 of ground state 
atom pairs calculated with the eigenfunction-based random 
phase approach (using 200 realizations for each J). 



Plotting the overlap, C^ n j versus eigenenergy, E s n J; of 

the electronic ground state Hamiltonian, H , yields os- 
cillations around a straight line, see Fig. [5] The slope 
multiplied by the Boltzmann constant represents the in- 
verse temperature, Convergence is defined by ap- 
proaching this line. A perfect straight lines means that 
a temperature can reliably be extracted from our numer- 
ical representation of the thermal ensemble. Decreasing 
oscillations around a straight line indicate that the tem- 
perature estimated from the thermally averaged overlaps 
is close to the actual temperature of the experiment. 

For the grid based random phase approach, cf. sec- 
tion 1111 Al the number of realizations, N, required to 
reach convergence was found to be much larger than 
the number of grid points Nr = 1024. As the number 
of realizations, N, is increased, the oscillations around 
the straight line decrease in Fig. [5] Note that for the 
eigenfunction-based random phase approach this con- 
sistency check with respect to temperature will give a 
straight line by construction. The free propagation ap- 
proach will approximate tightly the Boltzmann weighting 
of the scattering states. 

The initial thermal density of atom pairs, cf. Eq. 
is shown as a function of interatomic distance in Fig. [3] 
for the eigenfunction based random phase approach. For 
photoassociation, distances smaller than ~ 15 an are rele- 
vant. The thermal density is converged in this region by 
including rotational quantum numbers up to J = 200. 
The contribution of higher partial waves only ensures a 
constant density at large interatomic distances, see in- 
sert of Fig. [3] The long-distance part naturally converges 
very slowly but this is irrelevant for the dynamical cal- 
culations. The peak at short interatomic separations is 
due to bound levels, shape resonances and the classical 
turning point of the scattering states at the repulsive bar- 
rier of the potential: The difference between the black 
and yellow curves in Fig. [3] indicates the contribution of 




Grid based random phase 

Eigenfunction based random phase 
Excluding bound states 



Free propagation approach 
Eigenfunction based 
random phase 

Shape Resonance 



X L = 840 nm" 



50 100 150 

Rotational quantum number J 



FIG. 4. Population transfered to 1 Tl g state by the pump 
pulse, for different initial partial waves J, averaged over ran- 
dom phase wavefunctions (number of realizations: 200 for 
the 200 lowest eigenfunctions, 1030 in the grid-based random 
phase approach and 17 in the free propagation approach). 
The shaded part shows the contribution of bound levels in 
the eigenfunction based approach. The orange shaded part 
shows the contributions of the shape resonances. The statis- 
tical errors in these curves are estimated as a /mean ~ 0.11 
for the eigenvalue approach, a /mean ~ 0.09 for the grid- 
based approach and a /mean ~ 0.26 for the free propagation 
approach. 



bound levels, the difference between the yellow and the 
purple curve that of shape resonances. 

The population transferred from the initial incoherent 
ensemble to the 1 H g state by the pump pulse is shown in 
Figure |4j comparing the grid-based, eigenfunction-based 
and free propagation based random phase approach. The 
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thermal averaging procedure has been repeated for in- 
creasing initial rotational quantum number, J. That 
is, for each rotational barrier, grid-based, eigenfunction- 
based and Gaussian random phase wave functions are 
propagated in real time with the full, time-dependent 

Hamiltonian, H pump . Expectation values, such as the 
population of the 1 Ii g state after the pump pulse is 
over, are calculated for each random phase realization, 
k, or propagation time, t prop , respectively, and averaged 
over, including the rotational degeneracy factor J + 1, 
cf. Eqs. (H7J), (HU), and (gO]). The three approaches qual- 
itatively give the same result with a steep rise at low 
J-values, a peak at J = 55 and an exponential tail for 
J > 120. Each random phase approach represents a sta- 
tistical sampling of the photoassociation yield. For small 
iV the free propagation approach converged extremely 
fast. For refined results more sampling is required. The 
deviation of an expectation value from its mean scales 
as l/v^/V where N is the number of realizations. This 
was checked for J = 55 and J = 100 the pre-factor 
a /mean = s/ y^N) is estimated as s ~ 0.37 for the free 
propagation method s ~ 0.17 for the eigenvalue method 
and s ~ 0.30 for the grid based method. This makes the 
eigenvalue method converge fastest for large N. 

For low partial waves, the grid based (red curve) and 
the eigenfunction based (blue curve) random phase ap- 
proaches agree very well with each other. Once the bound 
levels are removed from the eigenfunction based approach 
(dashed purple curve), the eigenfunction-based approach 
agrees with the free propagation approach (green curve). 
This comparison allows for estimating the contribution 
of the bound levels. For J > 50 the ground state poten- 
tial does not support any bound levels due to the high 
centrifugal barrier. The total contribution of the bound 
part of the spectrum to the excitation of U g population 
amounts to about 15%. The free propagation approach 
does not account for the bound part since it builds a 
Boltzmann ensemble in the asymptotic, flat part of the 
potential, based on kinetic energy only. Relatively large 
deviations between the free propagation approach and 
the other methods are also observed for 50 < J < 65. 
This is partially attributed to shape resonances residing 
behind the centrifugal barrier. Due to the finite propaga- 
tion time (tp r0 p = 150 ps) in the free propagation method, 
only those resonances are captured that have a lifetime 
that is shorter than t prop . The contribution of long-lived 
shape resonances is absent since these resonances are not 
sampled by the free propagation approach. When the 
contributions of the shape resonances are added to the 
free propagation method (cyan curve) some of the dis- 
crepancy is removed. The differences between J = 60 to 
J = 90 are attributed to insufficient sampling of the free 
propagation method (only 17 realizations). At large J, 
the grid based random phase approach overestimates the 
exponential tail. This can be rationalized in terms of a 
biased sampling with too many points dedicated to high 
energies on the repulsive part of the centrifugal barrier. 

The role of shape resonances is further analyzed in 
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FIG. 5. Top panel: Population photoassociated out of shape 
resonances vs initial partial wave J. Bottom panel: Position 
in Kelvin (empty symbols) and lifetimes in nanoseconds (filled 
symbols) of the shape resonances vs partial wave J. 



Fig.[5]which displays the total photoassociation yield that 
originates from shape resonances (top panel) as well as 
the positions and lifetimes of the shape resonances (bot- 
tom panel). The shape resonances were calculated us- 
ing a complex absorbing potential^ Diagonalization of 
the non-Hermitian Hamiltonian was followed by propa- 
gation under the ground state Hamiltonian to get suf- 
ficiently accurate lifetimes. Shape resonances are found 
for J = 49, . . . , 62. For J = 49, 50, four shape resonances 
reside behind the centrifugal barrier, for J = 51, . . . , 55 
three, for J = 56, 57 two and for J = 58, . . . , 62 one. 
The positions of the resonances lie between 15 K and 
300 K and scale linearly with the rotational quantum 
number J (left-hand y-axis of the bottom panel of Fig. [5]). 
The lifetimes scale step-wise exponentially with J (and 
thus with energy) with the lowest-lying resonances hav- 
ing the longest lifetimes. Exponential scaling is expected 
if the lifetime is caused by tunneling out of the rotational 
barrier alone. The multi-exponential scaling can be at- 
tributed to the appearance of further resonances with the 
same rotational quantum number. For J = 49 ... , 57, 
i.e., for those partial waves where very long-lived reso- 
nances are observed, the contribution of the shape reso- 
nances to the photoassociation yield is fairly significant. 
It amounts to 23% for J = 50 and 9% for J = 55. This 
is not surprising since shape resonances represent quasi- 
bound states that are ideally suited for photoassociation. 
They give structure to the continuum of scattering states 
which otherwise is completely fiat at high temperature. 
This can be utilized for generation of coherence and con- 
trol. 

The missing contribution of photoassociation out of 
shape resonances in the free propagation approach ex- 
plains the discrepancy with the eigenfunction based 
method for the photoassociation yield in Fig. @] If the 
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FIG. 6. Top panel: Excited state purity, V R {tji na i), and co- 
herence measure, C e (t final), cf. Eqs. (|24|l and (|25|l . vs pulse 
intensity for the subensemble of photoassociated molecules. 
Bottom panel: Photoassociation yield for J = 55 as a func- 
tion of peak pulse intensity. 



contribution of shape resonances to the photoassociation 
yield is subtracted from the result of the eigenfunction 
based approach, the remaining error of 3% for J = 50 
and 10% for J = 55 is easily rationalized in terms of the 
statistical sampling. 

In conclusion, the free propagation approach is numer- 
ically most efficient method for a rough estimate. When 
further refinement is required the eigenvalue approach 
converges faster by a factor of 2. The free propagation 
approach excludes the bound part of the spectrum and 
the resonances, it comes with error bars of about 15%. 
If more accurate results are desired, the eigenfunction 
based random phase approach is the method of choice. 
The eigenfunction based method is also best suited to 
capture the contribution of bound states and shape res- 
onances to the photoassociation yield. 



VI. INTERPLAY OF COHERENCE AND 
PHOTOASSOCIATION YIELD 

In order to rationalize and quantify the generation 
of the coherence, that is observed experimentall y 28 ' 29 , 
out of a completely incoherent initial ensemble, we con- 
sider the change in quantum purity, Tr[p 2 ], due to 
the femtosecond laser pulses. The initial purity, deter- 
mined by the temperature, T = 1000 K, and density, 
p = 4.8 • 10 16 atoms/cm 3 , of the experiment, is estimated 
in terms of the purity of a single atom pair represented 
in our computation and the probability, p2, of finding 
two atoms in our computation volume, V. Given p and 
V = 4/37rr 3 = 6.33 • 10~ 21 cm 3 , we find p 2 = 9.2 • 10~ 8 . 
The purity of a single atom pair in our computation box 
is given by the vibrational purity weighted by the proba- 
bility, Pj, for occupation of the sub-ensemble with angu- 



lar momentum J. For an atom pair in their electronic 
ground state, the purity becomes V 9 = P%J2 1 p j^j 

where V 9 = Tr[(e~^ ) 2 }/Z 2 j, and and Zj are the 
Hamiltonian for vibrational motion and partition func- 
tion, respectively, of partial wave J. Evaluating Vj with 
freely propagated thermal random-phase wavefunctions, 
we obtain V 9 = 2.6 ■ 10~ 20 for the initial purity. While 
the effect of bound levels and shape resonances is not 
accounted for by the free propagation approach, their 
contribution to the initial purity is negligible. This is 
due to the fact that random phase wavefunctions for par- 
tial waves up to J = 800 need to be considered to ob- 
tain converged results. The bound levels and shape res- 
onances occur for small bands of J. Their weight within 
the large ensemble consisting of 800 partial waves is cor- 
respondingly small. For the ensemble of molecules in 
the electronically excited 1 I1 S state, the density opera- 
tor is given by p e = ^ZjPjp"j- Here, Pj is the prob- 
ability for occupation of the excited state sub-ensemble 
with angular momentum, J, cf. Fig. [5] p e j is the nor- 
malized density operator of the excited state J sub- 
ensemble which is constructed by incoherently averag- 
ing dyadic products of the propagated thermal wavefunc- 
tions, \ipn,j(tfi na i))(ip n ,j(tfi na i)\, over all random phases 
n and partial waves J (tfi na i denotes the time when the 
pump pulse is over). The purity for the excited state sub- 
ensemble is then given by V e = J2j( p j) 2 Tr [(Pj) 2 ]- We 
obtain a purity V e = 3.04 ■ 10~ 4 for the molecular sub- 
ensemble in the 1 H g excited state for the experimental 
pulse parameters. We thus observe a dramatic increase in 
the quantum purity, Tr[/3 2 ], induced by the femtosecond 
laser pulse. The underlying physical mechanism can be 
viewed as " Franck-Condon filtering" : for a given initial 
J value there is only a limited range of collision energies 
that allow the colliding pair to reach the Franck-Condon 
window for PA located at short internuclear distances-^. 

In order to obtain a quantitative estimate of the de- 
gree of distillation achieved by the fs-PA process, we 
have calculated the purity of ensemble of photoassoci- 
ated molecules in the ^Tlg state for a range of laser 
intensities, cf. Fig. \6\ For weak fields, the purity is 
roughly constant as a function of intensity and almost 
three times larger than the purity obtained for the in- 
tensity of 5 x 10 12 W/cm 2 used in the experiment. As 
intensity is increased, a sudden drop in the purity is 
observed which levels off at large intensities. We at- 
tribute this drop to power broadening for strong fields, 
which brings more atom pairs into the Franck-Condon 
window for PA. The purity increase with respect to the 
initial ensemble is nevertheless many orders of magnitude 
for both weak and strong field. To further analyze the 
generation of quantum coherence, it is possible to sepa- 
rate static and dynamic contributions, p = p stat + Pd yn - 
This is most easily achieved in the energy representation 
where the static (dynamic) part corresponds to the di- 
agonal (off-diagonal) matrix elements. The dynamical 
contributions are quantified by the coherence measure 
C = Tr[/3 2 , yn (i)]i£. Tracing the coherence measure of the 
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excited state, C e , as a function of laser intensity (red di- 
amonds in Fig. IH]) shows that most of the purity comes 
from the dynamical contribution (C e = 2.86 • 10~ 4 for the 
experimental intensity of 5 x 10 12 W/cm 2 ). 



VII. CONCLUSIONS 

We have described two-photon femtosecond photoas- 
sociation of magnesium atoms from first principles using 
state of the art ab initio methods and quantum dynami- 
cal calculations. Highly accurate potential energy curves 
were obtained using the coupled cluster method and a 
large basis set. Two-photon couplings and dynamic Stark 
shifts are important to correctly model the interaction of 
the atom pairs with the strong field of a femtosecond laser 
pulse. They were calculated within the framework of 
the equation of motion (linear response) coupled cluster 
method. The photoassociation dynamics were obtained 
by solving the Schrodingcr equation for all relevant par- 
tial waves, accounting for the laser-matter interaction in 
a non-perturbative way. 

An efficient numerical method has been developed to 
describe the incoherent thermal ensemble that is the ini- 
tial state for photoassociation at high temperatures. It is 
based on random phase wave packets which can be built 
from eigenfunctions of the grid, the Hamiltonian, or the 
kinetic energy. The latter can supply a rough estimate 
with ~ 17 realizations. This approach yields qualitatively 
correct results. It neglects, however, the contribution 
from bound levels and long-lived shape resonances and 
therefore comes with error bars of about 15%. The best 
compromise between high accuracy and convergence is 
found for the eigenfunction-based method where random 
phase realizations are built from the eigenfunctions of the 
electronic ground state Hamiltonian. About 200 partial 
waves and 200 realizations for each partial wave are re- 
quired for converged photoassociation dynamics. Time- 
dependent thermal averages are obtained by propagating 
each of the random phase basisfunctions and incoherently 
summing up all single expectation values. 

The random phase approach allows for constructing 
the thermal atom pair density as a function of inter- 
atomic separation for high temperatures. This is im- 
portant to highlight the difference between hot and cold 
photoassociatio n 42 ' 45 ' 50 . In the cold regime, the largest 
density is defined by the quantum reflection and resides 
in the long-distance, downhill part of the potential. The 
opposite is true in the hot regime: Here, the largest den- 
sity is found in the repulsive part of the ground state 
potential. This is due to the many partial waves that are 
thermally populated and the colliding atom pairs having 
sufficiently high kinetic energy to overcome the rotational 
barriers. For specific partial waves, shape resonances are 
found to be important. This is not surprising since they 
represent quasi-bound states that are ideally suited for 
photoassociation. At very low temperatures, most par- 
tial waves are frozen out and the scattering is almost 



exclusively s-wave. The role of the rotational quantum 
numbers J is less important in the electronically excited 
state but it is still detectable in form of quantum beats^ 
Both hot and cold photoassociation come with advan- 
tages as well as drawbacks. In the hot regime, molecules 
with much shorter bond length than in the cold regime 
are formed. However, the quantum purity and coherence 
of the created molecules is much larger in the cold regime 
where dynamical correlations exist prior to photoassoci- 
ation. These correlations indicate pre-entanglcment of 
the atom pair. Making a molecule corresponds to entan- 
gling two atoms, and photoassociation amounts to filter- 
ing out and entangled subensemble both in the hot and 
cold regime. 

Our study has opened up the possibility to study 
femotsecond photoassociation and its control at high 
temperatures and to investigate systematically the gen- 
eration of coherence out of an incoherent initial state. 
Future work will address the efficient theoretical descrip- 
tion of the probe step. 
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